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Abstract. - We investigate shock waves in the unitary Fermi gas by using the zero-temperature 
equations of superfluid hydrodynamics. We obtain analytical solutions for the dynamics of a 
localized perturbation of the uniform gas. These supersonic bright and subsonic dark solutions 
produce, after a transient time, an extremely large (divergent) density gradient: the shock. We 
calculate the time of formation of the shock and also simulate the space-time behavior of the waves 
by solving generalized hydrodynamic equations, which include a reliable dispersive regularization 
of the shock. We find that the shock spreads into wave ripples whose properties crucially depend 
on the chosen initial configuration. 



One of the basic problems in physics is how density per- 
turbations propagate through a material [l][2] . In addition 
to the well-known sound waves, there are shock waves 
characterized by an abrupt change in the density of the 
medium 1,2 . Shock waves are ubiquitous and have been 
studied in many different physical systems [H2] . Ten years 
ago shock waves have been experimentally observed also 
in atomic Bose- Einstein condensates (BECs) [2H2], and 
theoretically investigated in various BEC configurations 
[7HH] . Very recently the observation of nonlinear hydro- 
dynamic waves has been reported in the collision between 
two strongly interacting Fermi gas clouds of 6 Li atoms [TS] . 
The experiment shows the formation of density gradients, 
which are nicely reproduced by hydrodynamic equations 
with a phenomenological viscous term [15]. Nevertheless, 
the role of dissipation is questionable [16] since the ultra- 
cold unitary Fermi gas is noted as an example of an almost 
perfect fluid [17]. Indeed in Ref. [16] it has been shown, 
by solving zero-temperature time-dependent Bogolibov-de 
Gennes equations, that the viscous term is not necessary 
to reproduce the experimental results of Ref. [IS] . 

Here we investigate the formation and dynamics of 
shock waves in the dilute and ultracold unitary (diver- 
gent inter-atomic scattering length) Fermi gas by using 
the zero-temperature equations of superfluid hydrodynam- 
ics pp. At zero temperature fermionic superfluids in the 
BCS-BEC crossover can be modelled by hydrodynamic 
equations [TH] and their generalizations with a gradient 
term |19) that induces a dispersive regularization of the 
shock. In this paper we obtain analytical solutions for 



the dynamics of a localized perturbation of the uniform 
gas. We calculate the supersonic (or subsonic) velocity of 
propagation of these bright (or dark) perturbations. We 
show that bright perturbations evolve towards a shock- 
wave front, while dark perturbations produce the shock in 
their back, and we calculate the period T s of formation of 
the shock. In addition, we study the space-time behavior 
of the shock waves beyond this characteristic time T s by 
including a reliable quantum correction in the hydrody- 
namic equations [15]. In fact, according to the two-fluid 
model of Landau [T], the viscous term acts only on the 
normal component of the fluid, and at zero temperature 
the normal component is zero. Morover, recent theoretical 
microscopic calculations [20] suggest that the viscosity of 
the unitary Fermi gas is extremely small at very low tem- 
peratures because the transverse current does not couple 
to collective modes. By solving numerically these general- 
ized equations of the unitary Fermi gas we show that the 
shock-wave front spreads into wave ripples whose proper- 
ties crucially depend on the "brightness" (bright or dark) 
of the chosen initial configuration. We expect that our re- 
sults are reliable when the normal density (with its viscous 
term) is quite small, i.e. for a temperature T much smaller 
than the critical temperature T c of the normal-superfluid 
transition. For the unitary Fermi gas T c ~ 0.2 Tp, where 
Tf = £f /kB is the Fermi temperature with ks the Boltz- 
mann constant, ep = ?i 2 (37r 2 n) 2 / 3 /(2m) the bulk Fermi 
energy of the ideal Fermi gas, n the total density, and m 
the mass of each atomic fermion. According to Ref. [21] 
for T < 0.1 Tp the normal fraction is below 10% and, 
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in practice, in this range of temperatures our approach is 
fully justified. 

At zero temperature the low-energy collective dynam- 
ics of a Fermi superfluid of neutral and dilute atoms at 
unitarity can be described by the equations of irrotational 
and inviscid hydrodynamics 



superfluid, given by 
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where n(r,t) is the total density of the superfluid, v(r, t) 
is its velocity field, and 



'2m' 



(3) 



is the bulk chemical potential of the system, with £ ~ 
0.4 a universal parameter |18j . Here we are supposing a 
balanced system, namely the same number of fermions in 
the two components of the spin (a \). The term U(r) 
models the external potential which traps the atoms. 

Let us consider the unitary Fermi gas with constant den- 
sity n with U(r) = 0. Experimentally this configuration 
can be obtained with a very large square-well potential 
(or a similar external trapping), such that in the model 
one can effectively impose periodic boundary conditions 
instead of the vanishing ones. A density variation along 
the z axis with respect to the uniform configuration n can 
be experimentally created by using a blue-detuned (bright 
perturbation) or a red-detuned (dark perturbation) laser 
beam |22j . In practice, we perform the following factor- 
ization 

n(r) = nx(x,y)m(z) , (4) 



by imposing also 



such that 



n±(x,y) = n± 
n \\( z ) = n\\ p( z ) 

n(r) — n p(z) 
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where n = n±_fi\\ , and p(z, t) is the relative density, i.e. the 
localized axial modification with respect to the uniform 
density n. We impose periodic boundary conditions along 
the z axis, namely p(z — L z ,t) = p(z = —L z ,t), with 2L Z 
the axial-domain length. We set v(r, t) = (0, 0, v(z, t)) 
with v(z,t) the velocity field such that v(z — L z ,t) = 
v(z = —L z ,t). Moreover we impose that the initial 
localized wave packet satisfies the boundary conditions 
p(z = ±L z ,t = 0) = 1 and v(z = ±L z ,t = 0) = 0. 
Because the dimensional reduction is done assuming the 
uniformess in x,y directions, we shall consider the propa- 
gation of a plane wave along the z axis. 

Inserting Eq. (0 into Eqs. (fTJ and ([2]) one finds the 
ID hydrodynamic equations for the axial dynamics of the 



p + vp + v p ~ , 



(8) 
(9) 



where dots denote time derivatives, primes space deriva- 
tives, and 



ci s (p) =c s p 



1/3 



(10) 



is the local sound velocity, with c s = q s (1) = \/£,/3vf the 
bulk sound velocity, v F — y / 2ep/m is bulk Fermi velocity 

and 6f = ^(3tt 2 «) 2 ^ 3 the bulk Fermi energy. 

The bulk sound velocity c s is the speed of propagation 
of a small perturbation p(z, t) with respect to the uniform 
superfluid of density n. In fact, setting p{z, t) = l+p(z,t), 
with p(z,t) <C 1 and v(z,t) of the same order of p(z,t), 
from the linearization of Eqs. ((5J and ([9]) we get the fa- 
miliar linear wave equation 



dt 2 
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for p{z, t) and a similar equation for v{z, t). Modelling the 
initial perturbation with a Gaussian shape, i.e. 



p{z,{)) = 2r 1 e - z ' 2 l^ 



(12) 



one finds [TJ[2] from the linearized equations 

p( Z , t) = n e -(^stf/(2**) + ^ e -( Z+ c s tf/(2**) ; (13) 

with initial condition p(z,t = 0) = 0. Thus, for the con- 
servation of the linear momentum, the initial wave packet 
splits into two waves travelling in opposite directions with 
the speed of sound c s . Obviously Eq. ([rS]) is reliable only 
if \r]\ <^ 1. As expected, a small (infinitesimal) perturba- 
tion gives rise to sound waves. 

We can find wave solutions of Eqs. ([8]) and (|9|) with 
a generic initial density profile by following the approach 
described by Landau and Lifshitz pQ. By supposing that 
the velocity v depends explicitly on the density p, i.e. v — 
v(p(z,t)), one has v = j^p, v' = j^p' ■ We now impose 
that the two hydrodynamic equations reduce to the same 
hyperbolic equation 



p + c{p)p' = , 



where 



c{p) = v{p) 
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It is quite easy to verify that, given a initial condition 
F(z) — p(z,t — 0) for the density profile, the time- 
dependent solution p(z,t) of the hyperbolic equation (JTJJ 
satisfies the following implicit, but algebraic, equation: 



p(z,t) = F{z - c(p(z,t))t) 



(16) 
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Fig. 1: Properties of the shock waves. Upper panel: Mach 
number M as a function of the amplitude r\ of the pertur- 
bation (solid line). For M < 1 there are subsonic and dark 
( — 1/2 < 77 < 0) waves, while for M > 1 there are supersonic 
and bright (77 > 0) waves. The dotted line separates the two 
regimes. Lower panel: period T s of formation (breaking time) 
of the shock-wave front as a function of the amplitude 77 of the 
perturbation. T a is in units of <r/c s , where a is the width of 
the perturbation and c s — -J £/3vf is the bulk speed of sound, 
with vf the Fermi velocity. 



To determine the local velocity of propagation c(p), 
which is not equal to the sound velocity ci s (p), we observe 
that from Eqs. (fTS]) we get 



dv _ c is (p) 
dp p 



(17) 



After separation of variables, and imposing that at infinity 
the density is equal to one and the velocity field is zero, 
we finally get 



v{p) = ±3c s (p 1 / 3 - 1) 



(18) 



The velocity c(p) follows directly from the velocity v{p) 
by using Eqs. (IT51) . One finds c(p) — v(p)±ci s (p), namely 



c(p) = ±c s (V /3 ~ 3) 



(19) 



In conclusion, we have found that the density field p(z,t) 
satisfies the implicit algebraic equation (|16[) with c{p) 
given by Eq. (fT§|) . Note that a similar result, but with a 
very different local velocity, has been obtained by Damski 
[8 for the ID BEC. 

Eqs. (fTB) and (fT9|) contain the dynamics of the two 
waves propagating to the left and to the right with ini- 
tial condition (|12j) . Some properties characterizing the 
dynamics can be extracted from these equations. First of 
all the two travelling waves have symmetric shapes during 
the time evolution. In addition, both amplitude and ve- 
locity of the extrema (maxima or minima, depending on 



the sign of 77) of the two waves are practically constant 
during time evolution. In particular, the amplitude of the 
extrema is given by A(rj) = 1 + 77 while the velocity of the 
extrema reads 



V{rf) = c(l + r/) = ±c s (4(1 + 7/) 1 / 3 - 3) 



(20) 



Notice that taking 77 = the velocity of the impulse ex- 
trema reduces to the sound velocity: V(0) = c(l) = c s = 
y/^/Svp- Moreover, bright perturbations (77 > 0) move 
faster than dark ones (77 < 0), and the Mach number 
M — V(t))/V(0) of these perturbations in the unitary 
Fermi gas is simply 



M = 4(1 + 7/) 1/3 - 3 . 



(21) 



For M > 1, which means 7/ > (bright perturbation), one 
has supersonic waves, while for < M < 1, which means 
77 < (dark perturbation), one has subsonic waves. In 
the upper panel of Fig. Q] we plot the Mach number M 
as function of the amplitude 77 of the perturbation. Note 
that since 2?7 is the amplitude of the initial condition, see 
Eq. (|T2"|) . the region 77 < —1/2 is unphysical. 

Let us consider a bright perturbation (77 > 0) moving to 
the right. The speed of impulse maximum V(j]) is bigger 
than the speed of its tails V(0). As a result the impulse 
self-steepens in the direction of propagation and a shock 
wave front takes place. The breaking-time T s required 
for such a process can be estimated as follows: the shock 
wave front appears when the distance difference traveled 
by lower and upper impulse parts is equal to the impulse 
half-width 2a, namely [V(rj) - V(0)}T s = 2a. This, by 
using Eqs. (fl9"]l and Eq. (fT0|). gives 



T, = 



2c s ((l + 7?) 1 /3_l) 



(22) 



In the case of a dark perturbation (77 < 0) the tails of 
the wave packet move faster than the impulse minimum. 
The shock appears in the back of the travelling wave, and 
the period of shock formation is simply T s = 2<r/(V(0) — 
V{rj)). In the lower panel of Fig. Q]we plot the period T s 
as a function of the amplitude 77 of the perturbation. The 
figure shows that as 77 goes to zero the period T s goes to 
infinity; in fact, in this limit the shock wave reduces to a 
sonic wave (sound wave) which does not produce a shock. 

After the formation of the shock Eqs. ([1]) and @ are 
not reliable because their exact solutions given of Eqs. 
(fT6| and (fl"9)) are no more single-valued. To overcome 
this difficulty we include a gradient quantum term in the 
hydrodynamic equations, which become 
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^n + V-(nv) = 0, (23) 
+ p(n) + f/(r)l = , (24) 



We stress that at zero temperature the simplest regular- 
ization process of the shock is a purely dispersive quantum 
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Fig. 2: Time evolution of supersonic shock waves. Initial 
condition with ajlp — 18 and rj = 0.3. The curves give 
the relative density profile p(z) at subsequent frames, where 

If ~ \J h 2 1 (meF) is the Fermi length and ljf = (.f /Ti is the 
Fermi frequency. 



gradient term. Historically, the gradient term with A was 
introduced by von Weizsacker to treat surface effects in 
nuclei |23) . This approach has been adopted for quantum 
hydrodynamics of electrons by March and Tosi [24], and 
also by Zaremba and Tso [35]. In the study of the BCS- 
BEC crossover the gradient term has been considered by 
various authors [53] ■ The choice of the parameter A in Eqs. 
(123)) and (|2~Ij) is still debated, here we choose A = 1/4, a 
value which gives good agreement with Monte Carlo cal- 
culations at zero and finite temperature (for details see 
[19 1 121 ] ). Moreover we set £ = 0.4. 

We expect that Eqs. ([23]) and (|24]l are reliable to study 
the long-time dynamics of shock waves in the ultracold 
unitary Fermi gas. It is well known that, according to the 
two-fluid model and the Landau's criterion of superfluidity 
PP, above a critical temperature v c a normal component 
with a dissipative term appears in the fluid [1 . As dis- 
cussed in [57], for the unitary Fermi gas one has v c ~ c s . 
Nevertheless, at ultrcold temperature the normal compo- 
nent is negligible and also the shear viscosity [17j[20]. For 
these reasons at zero temperature the shock waves are dis- 
persive and not dissipative [16) . 

Eqs. (123)) and (|2"^1) can be formally written (for any 
value of A, also A = 0) in terms of a Galilei-invariant 
nonlinear Schrodinger equation [115]. Setting U(r) = 
and using Eqs. ([3]) and ([7]) we easily get from Eqs. ([2"3"]) 
and (|24]) a ID nonlinear Schodinger equation. We solve 
this equation by using a Crank-Nicolson finite-difference 
predictor-corrector algorithm [2 8) with the initial condi- 
tion given by Eq. (fl"2"j) and v(z,t = 0) = 0. In fact, as 
also shown by Damski [5], we have verified that the ini- 
tial velocity field v(p(z,t = 0)) and v(z,t = 0) = give 
practically the same time evolution. 



In Fig. [2] we plot the time evolution of supersonic 
shock waves obtained with cr/lp = 18 and rj = 0.3, with 

J h 2 1 (mep) the Fermi length of the bulk system. The fig- 
ure displays the density profile p{z) at subsequent times. 
Note the splitting on the initial bright wave packet into 
two bright travelling waves moving in opposite directions. 
As previously discussed, there is a deformation of the two 
waves with the formation of a quasi-horizontal shock- wave 
front. Eventually, this front spreads into wave ripples. 
There is no qualitative difference with respect to a Bose- 
Einstein condensate [8] in the physical manifestation of 
supersonic shock waves in the zero-temperature unitary 
Fermi gas. Nevertheless, due to the very different equa- 
tion of state, there are large quantitative differences. Our 
numerical analysis confirms that the breaking time T s de- 
creases by increasing the amplitude 77, while the velocity V 
of the maxima of the travelling waves increases by increas- 
ing 77. There is a good agreement between our analytical 
formulas, Eqs. (|20|) and (f22|) . and simulations: the rela- 
tive difference is within 5% for the velocity V and within 
20% for the breaking time T s . 

In Fig. [3] we plot the time evolution of subsonic shock 
waves obtained again with <j/If — 18 and r\ — —0.2. Also 
in this case the figure shows the splitting on the initial 
dark wave packet into two dark travelling waves moving 
in opposite directions. But here, as expected, the quasi- 
horizontal shock appears in the back side of the travelling 
waves. Our simulations show that for 77 < the wave rip- 
ples which appear at the breaking time T s are always dark, 
i.e. they never exceed the bulk density (compare wave rip- 
ples of Fig. [5] with those of Fig. [3]). Note that dark shock 
waves have been studied long ago [29] in a different phys- 
ical context: the discrete nonlinear Schrodinger equation. 
In that case the wave ripples can exceed the bulk density, 
probably due to the discrete nature of the Schrodinger 
equation. Also for dark shock waves our analytical pre- 
dictions on velocity V of the minima and breaking time 
T s are quite accurate with respect to numerical findings 
(similar relative differences of runs with rj > 0). 

In conclusion, we have shown that at very low temper- 
atures the unitary Fermi gas admits supersonic and sub- 
sonic shock waves, for which we have developed analytical 
and numerical results. Our predictions suggest a much 
cleaner method to produce shock waves with respect to 
the recent experiment [15] based on the collision of two 
6 Li atomic clouds. The shape of these waves changes dur- 
ing the time evolution giving rise to a shock-wave front 
at a characteristic breaking time. We have determined 
the Mach number of these travelling waves as a function 
of the perturbation amplitude, showing that supersonic 
bright and subsonic dark waves behave quite differently. 

The author thanks Sadahn Kumar Adhikari, Tilman 
Enss, Boris Malomed, Enzo Orlandini, Mario Salerno, and 
Flavio Toigo for useful suggestions and discussions. 
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Fig. 3: Time evolution of subsonic shock waves. Initial con- 
dition with a /If = 18 and n = —0.2. The curves give 
the relative density profile p(z) at subsequent frames, where 

If — \J ti 2 1 (mei?) is the Fermi length and ljf = £F/fi is the 
Fermi frequency. 
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